###############################################################
###############################################################
###                                                         ###
### Replication Code for Krogslund, Choi, & Poertner (2014) ###
###                                                         ###
###############################################################
###############################################################

### Set main working directory
main="/Users/ChristopherKrogslund/Documents/Papers/Krogslund, Choi, Poertner (2014)"
setwd(main)

### Clear workspace
rm(list=ls())

### Load required packages and functions 
library("ggplot2")
library("plyr")
library("reshape2")
library("QCA")
library("multicore")
source(file="R Functions/fsQCA_sim_inclcut.R")
source(file="R Functions/fsQCA_sim.R")
source(file="R Functions/calibrate_samford.R")
source(file="R Functions/fsQCA_anchor_samford.R")
source(file="R Functions/calibrate_ahn_lee.R")
source(file="R Functions/fsQCA_anchor_ahn_lee.R")
source(file="R Functions/fsQCA_anchor_ka.R")
source(file="R Functions/fsQCA_random.R")


#####################################################
### Replication & Simulation for Ahn & Lee (2012) ###
#####################################################

### Set working sub-directory
setwd("./Ahn & Lee (2012)")

### Import the data set
data <- read.csv(file="CWS_data_extended_KOR.csv")

### Extract Korea data
kor.data=data[which(data$ID=="KOR"&data$YEAR>=1964),]

### Create WELGDP (taken to be SSEXP over GDPNCU, or social exp as % of current GDP in mil NCU)
kor.data["WELGDP"]=kor.data$SSEXP/kor.data$GDPNCU/10

### Isolate the dependent variable
depvar=kor.data[, c("YEAR","WELGDP")]

### Create LTPRD (dummy for Kim Dae-Jung administration (leftgov))
kor.data["LTPRD"]=NA
for(i in 1:nrow(kor.data)){
  if(kor.data$YEAR[i]<1998){
    kor.data$LTPRD[i]=0
  } else(kor.data$LTPRD[i]=1)
}

### Create ELDERLY (percentage of the population 65+)
kor.data["ELDERLY"]=kor.data$PO65NEW/kor.data$POPNEW*100

### Create correct DINVOC 
kor.data["DINVOC"]=kor.data$DINVIC-kor.data$DINVOC

### Extract the independent variables
indepvar.1=kor.data[,c("YEAR", "CGDP", "ELDERLY", "STUNEMR", "LTPRD", "DINVOC")]

### Calibrate the raw data using given achor points
raw.data=join(x=depvar, y=indepvar.1, type="left")
WELGDP=c(1,3,7)
CGDP=c(2104, 9560, 20000)
ELDERLY=c(3.39,5.34, 9.5)
STUNEMPR=c(2.4, 3.75, 5)
DINVOC=c(-1.5, -0.02, 0.44)
user.thresholds=data.frame(WELGDP, CGDP, ELDERLY, STUNEMPR, DINVOC)
raw.data=raw.data[,-c(1)]
raw.data=na.exclude(raw.data)
repdata.cal=calibrate.ahn.lee(x=user.thresholds)

### Set working sub-directory
setwd("./Results")

### Simulate Model #1
sim.model.1=fsQCA.sim(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "DINVOC", "LTPRD"), min.freq=1, max.freq=10, reps=3000, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0)
write.table(x=sim.model.1$results, file="model_1_sim.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_sim.pdf", plot=sim.model.1$plot, scale=1.85)

### Simulate Model #2
sim.model.2=fsQCA.sim(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP"), min.freq=1, max.freq=10, reps=3000, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0)
write.table(x=sim.model.2$results, file="model_2_sim.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_sim.pdf", plot=sim.model.2$plot, scale=1.85)

### Simulate Model #3
sim.model.3=fsQCA.sim(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC"), min.freq=1, max.freq=10, reps=3000, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0)
write.table(x=sim.model.3$results, file="model_3_sim.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_sim.pdf", plot=sim.model.3$plot, scale=1.85)

### Simulate Model #4
sim.model.4=fsQCA.sim(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "LTPRD"), min.freq=1, max.freq=10, reps=3000, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0)
write.table(x=sim.model.4$results, file="model_4_sim.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_sim.pdf", plot=sim.model.4$plot, scale=1.85)

### Simulate Model #5
sim.model.5=fsQCA.sim(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC", "LTPRD"), min.freq=1, max.freq=10, reps=3000, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0)
write.table(x=sim.model.5$results, file="model_5_sim.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_sim.pdf", plot=sim.model.5$plot, scale=1.85)

### Simulate Model #6
sim.model.6=fsQCA.sim(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "LTPRD"), min.freq=1, max.freq=10, reps=3000, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0)
write.table(x=sim.model.6$results, file="model_6_sim.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_sim.pdf", plot=sim.model.6$plot, scale=1.85)


### Simulate Model #1 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.1.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "DINVOC", "LTPRD"), n.cut=1, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.1.error$results
write.table(x=results, file="model_1_ncut1_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut1_error.pdf", plot=sim.model.1.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_1_ncut1_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut1_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #2 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.2.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP"), n.cut=1, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.2.error$results
write.table(x=results, file="model_2_ncut1_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut1_error.pdf", plot=sim.model.2.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_2_ncut1_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut1_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #3 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.3.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC"), n.cut=1, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.3.error$results
write.table(x=results, file="model_3_ncut1_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut1_error.pdf", plot=sim.model.3.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_3_ncut1_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut1_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #4 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.4.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "LTPRD"), n.cut=1, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.4.error$results
write.table(x=results, file="model_4_ncut1_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut1_error.pdf", plot=sim.model.4.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_4_ncut1_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut1_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #5 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.5.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC", "LTPRD"), n.cut=1, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.5.error$results
write.table(x=results, file="model_5_ncut1_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut1_error.pdf", plot=sim.model.5.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_5_ncut1_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut1_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #6 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.6.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "LTPRD"), n.cut=1, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.6.error$results
write.table(x=results, file="model_6_ncut1_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut1_error.pdf", plot=sim.model.6.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_6_ncut1_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut1_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #1 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.1.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "DINVOC", "LTPRD"), n.cut=2, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.1.error$results
write.table(x=results, file="model_1_ncut2_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut2_error.pdf", plot=sim.model.1.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_1_ncut2_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut2_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #2 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.2.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP"), n.cut=2, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.2.error$results
write.table(x=results, file="model_2_ncut2_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut2_error.pdf", plot=sim.model.2.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_2_ncut2_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut2_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #3 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.3.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC"), n.cut=2, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.3.error$results
write.table(x=results, file="model_3_ncut2_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut2_error.pdf", plot=sim.model.3.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_3_ncut2_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut2_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #4 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.4.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "LTPRD"), n.cut=2, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.4.error$results
write.table(x=results, file="model_4_ncut2_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut2_error.pdf", plot=sim.model.4.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_4_ncut2_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut2_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #5 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.5.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC", "LTPRD"), n.cut=2, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.5.error$results
write.table(x=results, file="model_5_ncut2_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut2_error.pdf", plot=sim.model.5.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_5_ncut2_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut2_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #6 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.6.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "LTPRD"), n.cut=2, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.6.error$results
write.table(x=results, file="model_6_ncut2_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut2_error.pdf", plot=sim.model.6.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_6_ncut2_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut2_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #1 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.1.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "DINVOC", "LTPRD"), n.cut=3, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.1.error$results
write.table(x=results, file="model_1_ncut3_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut3_error.pdf", plot=sim.model.1.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_1_ncut3_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut3_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #2 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.2.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP"), n.cut=3, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.2.error$results
write.table(x=results, file="model_2_ncut3_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut3_error.pdf", plot=sim.model.2.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_2_ncut3_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut3_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #3 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.3.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC"), n.cut=3, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.3.error$results
write.table(x=results, file="model_3_ncut3_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut3_error.pdf", plot=sim.model.3.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_3_ncut3_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut3_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #4 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.4.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "LTPRD"), n.cut=3, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.4.error$results
write.table(x=results, file="model_4_ncut3_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut3_error.pdf", plot=sim.model.4.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_4_ncut3_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut3_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #5 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.5.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC", "LTPRD"), n.cut=3, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.5.error$results
write.table(x=results, file="model_5_ncut3_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut3_error.pdf", plot=sim.model.5.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_5_ncut3_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut3_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #6 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.6.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "LTPRD"), n.cut=3, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.6.error$results
write.table(x=results, file="model_6_ncut3_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut3_error.pdf", plot=sim.model.6.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_6_ncut3_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut3_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #1 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.1.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "DINVOC", "LTPRD"), n.cut=4, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.1.error$results
write.table(x=results, file="model_1_ncut4_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut4_error.pdf", plot=sim.model.1.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_1_ncut4_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut4_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #2 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.2.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP"), n.cut=4, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.2.error$results
write.table(x=results, file="model_2_ncut4_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut4_error.pdf", plot=sim.model.2.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_2_ncut4_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut4_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #3 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.3.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC"), n.cut=4, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.3.error$results
write.table(x=results, file="model_3_ncut4_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut4_error.pdf", plot=sim.model.3.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_3_ncut4_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut4_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #4 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.4.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "LTPRD"), n.cut=4, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.4.error$results
write.table(x=results, file="model_4_ncut4_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut4_error.pdf", plot=sim.model.4.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_4_ncut4_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut4_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #5 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.5.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC", "LTPRD"), n.cut=4, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.5.error$results
write.table(x=results, file="model_5_ncut4_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut4_error.pdf", plot=sim.model.5.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_5_ncut4_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut4_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #6 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.6.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "LTPRD"), n.cut=4, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.6.error$results
write.table(x=results, file="model_6_ncut4_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut4_error.pdf", plot=sim.model.6.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_6_ncut4_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut4_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #1 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.1.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "DINVOC", "LTPRD"), n.cut=5, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.1.error$results
write.table(x=results, file="model_1_ncut5_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut5_error.pdf", plot=sim.model.1.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_1_ncut5_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut5_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #2 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.2.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP"), n.cut=5, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.2.error$results
write.table(x=results, file="model_2_ncut5_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut5_error.pdf", plot=sim.model.2.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_2_ncut5_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut5_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #3 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.3.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC"), n.cut=5, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.3.error$results
write.table(x=results, file="model_3_ncut5_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut5_error.pdf", plot=sim.model.3.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_3_ncut5_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut5_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #4 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.4.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "LTPRD"), n.cut=5, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.4.error$results
write.table(x=results, file="model_4_ncut5_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut5_error.pdf", plot=sim.model.4.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_4_ncut5_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut5_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #5 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.5.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC", "LTPRD"), n.cut=5, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.5.error$results
write.table(x=results, file="model_5_ncut5_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut5_error.pdf", plot=sim.model.5.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_5_ncut5_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut5_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #6 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.6.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "LTPRD"), n.cut=5, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.6.error$results
write.table(x=results, file="model_6_ncut5_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut5_error.pdf", plot=sim.model.6.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_6_ncut5_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut5_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #1 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.1.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "DINVOC", "LTPRD"), n.cut=6, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.1.error$results
write.table(x=results, file="model_1_ncut6_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut6_error.pdf", plot=sim.model.1.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_1_ncut6_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut6_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #2 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.2.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP"), n.cut=6, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.2.error$results
write.table(x=results, file="model_2_ncut6_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut6_error.pdf", plot=sim.model.2.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_2_ncut6_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut6_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #3 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.3.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC"), n.cut=6, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.3.error$results
write.table(x=results, file="model_3_ncut6_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut6_error.pdf", plot=sim.model.3.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_3_ncut6_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_ncut6_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #4 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.4.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "LTPRD"), n.cut=6, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.4.error$results
write.table(x=results, file="model_4_ncut6_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut6_error.pdf", plot=sim.model.4.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_4_ncut6_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_ncut6_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #5 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.5.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC", "LTPRD"), n.cut=6, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.5.error$results
write.table(x=results, file="model_5_ncut6_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut6_error.pdf", plot=sim.model.5.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_5_ncut6_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_ncut6_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Simulate Model #6 with error in the cross-over point (0.00,0.01, 0.025, 0.05)
sim.model.6.error=fsQCA.anchor.ahn.lee(raw.data=raw.data, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "LTPRD"), n.cut=6, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.85, given.incl.cut1=0.85, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.6.error$results
write.table(x=results, file="model_6_ncut6_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut6_error.pdf", plot=sim.model.6.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_6_ncut6_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_ncut6_error_measure.pdf", plot=sensitivity.plot, scale=1.25)


### Simulate Model #1 with Random Variable
sim.model.1.random=fsQCA.random(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "DINVOC", "LTPRD"), min.freq=1, max.freq=10, reps=100, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, random.reps=100, min=0)
write.table(x=sim.model.1.random$results, file="model_1_sim_random.csv", sep=",", row.names=FALSE)
percent=c(sim.model.1.random$percent.pos, sim.model.1.random$percent.neg, sim.model.1.random$percent.all)
percent.label=c("Percent Positive", "Percent Negative", "Percent All")
percents=data.frame(percent.label, percent)
write.table(x=percents, file="model_1_sim_random_percents.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_sim_random.pdf", plot=sim.model.1.random$plot, scale=1.85)

### Simulate Model #2 with Random Variable
sim.model.2.random=fsQCA.random(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP"), min.freq=1, max.freq=10, reps=100, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, random.reps=100, min=0)
write.table(x=sim.model.2.random$results, file="model_2_sim_random.csv", sep=",", row.names=FALSE)
percent=c(sim.model.2.random$percent.pos, sim.model.2.random$percent.neg, sim.model.2.random$percent.all)
percent.label=c("Percent Positive", "Percent Negative", "Percent All")
percents=data.frame(percent.label, percent)
write.table(x=percents, file="model_2_sim_random_percents.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_sim_random.pdf", plot=sim.model.2.random$plot, scale=1.85)

### Simulate Model #3 with Random Variable
sim.model.3.random=fsQCA.random(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC"), min.freq=1, max.freq=10, reps=100, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, random.reps=100, min=0)
write.table(x=sim.model.3.random$results, file="model_3_sim_random.csv", sep=",", row.names=FALSE)
percent=c(sim.model.3.random$percent.pos, sim.model.3.random$percent.neg, sim.model.3.random$percent.all)
percent.label=c("Percent Positive", "Percent Negative", "Percent All")
percents=data.frame(percent.label, percent)
write.table(x=percents, file="model_3_sim_random_percents.csv", sep=",", row.names=FALSE)
ggsave(filename="model_3_sim_random.pdf", plot=sim.model.3.random$plot, scale=1.85)

### Simulate Model #4 with Random Variable
sim.model.4.random=fsQCA.random(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "LTPRD"), min.freq=1, max.freq=10, reps=100, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, random.reps=100, min=0)
write.table(x=sim.model.4.random$results, file="model_4_sim_random.csv", sep=",", row.names=FALSE)
percent=c(sim.model.4.random$percent.pos, sim.model.4.random$percent.neg, sim.model.4.random$percent.all)
percent.label=c("Percent Positive", "Percent Negative", "Percent All")
percents=data.frame(percent.label, percent)
write.table(x=percents, file="model_4_sim_random_percents.csv", sep=",", row.names=FALSE)
ggsave(filename="model_4_sim_random.pdf", plot=sim.model.4.random$plot, scale=1.85)

### Simulate Model #5 with Random Variable
sim.model.5.random=fsQCA.random(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "DINVOC", "LTPRD"), min.freq=1, max.freq=10, reps=100, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, random.reps=100, min=0)
write.table(x=sim.model.5.random$results, file="model_5_sim_random.csv", sep=",", row.names=FALSE)
percent=c(sim.model.5.random$percent.pos, sim.model.5.random$percent.neg, sim.model.5.random$percent.all)
percent.label=c("Percent Positive", "Percent Negative", "Percent All")
percents=data.frame(percent.label, percent)
write.table(x=percents, file="model_5_sim_random_percents.csv", sep=",", row.names=FALSE)
ggsave(filename="model_5_sim_random.pdf", plot=sim.model.5.random$plot, scale=1.85)

### Simulate Model #6 with Random Variable
sim.model.6.random=fsQCA.random(dataset=repdata.cal, depvar=c("WELGDP"), indepvars=c("ELDERLY", "STUNEMR", "CGDP", "LTPRD"), min.freq=1, max.freq=10, reps=100, explain=1, given.incl.cut0=0.85, given.incl.cut1=0.85, random.reps=100, min=0)
write.table(x=sim.model.6.random$results, file="model_6_sim_random.csv", sep=",", row.names=FALSE)
percent=c(sim.model.6.random$percent.pos, sim.model.6.random$percent.neg, sim.model.6.random$percent.all)
percent.label=c("Percent Positive", "Percent Negative", "Percent All")
percents=data.frame(percent.label, percent)
write.table(x=percents, file="model_6_sim_random_percents.csv", sep=",", row.names=FALSE)
ggsave(filename="model_6_sim_random.pdf", plot=sim.model.6.random$plot, scale=1.85)


###################################################
### Replication & Simulation for Samford (2010) ###
###################################################

### Set working sub-directory
setwd(main)
setwd("./Samford (2010)")

### Import the data set
data <- read.csv(file="Samford2010set.csv")

### Isolate dependent and independent variables used in the model
repdata.cal=data[,c(5:12)]

### Prepare raw data for error simulations
raw.data=read.csv(file="samford_raw.csv")
raplib=c(5,25,50)
manufac=c(10,20,30)
devalu=c(50,100,200)
hyperinf=c(10,50,100)
grostron=c(2,4,6)
groweak=c(0, -1.5, -3)
switcher=c(0, 0.5, 1)
execunco=c(6,4,2)
user.thresholds=data.frame(raplib, manufac, devalu, hyperinf, grostron, groweak, switcher, execunco)

### Set working sub-directory
setwd("./Results")

### Simulate Model
sim.model=fsQCA.sim(dataset=repdata.cal, depvar=c("raplib"), indepvars=c("manufac", "devalu", "hyperinf", "grostron", "groweak", "switcher", "execunco"), min.freq=1, max.freq=10, reps=3000, explain=1, given.incl.cut0=0.65, given.incl.cut1=0.65, min=0)
write.table(x=sim.model$results, file="model_sim.csv", sep=",", row.names=FALSE)
ggsave(filename="model_sim.pdf", plot=sim.model$plot, scale=1.85)


### Simulate the model with error in the cross-over point (0, 0.01, 0.025, 0.05)
sim.model.error=fsQCA.anchor.samford(raw.data=raw.data, depvar=c("raplib"), indepvars=c("manufac", "devalu", "hyperinf", "grostron", "groweak", "switcher", "execunco"), n.cut=1, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.65, given.incl.cut1=0.65, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.error$results
write.table(x=results, file="model_ncut1_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_ncut1_error.pdf", plot=sim.model.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_ncut1_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_ncut1_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

sim.model.error=fsQCA.anchor.samford(raw.data=raw.data, depvar=c("raplib"), indepvars=c("manufac", "devalu", "hyperinf", "grostron", "groweak", "switcher", "execunco"), n.cut=2, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.65, given.incl.cut1=0.65, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.error$results
write.table(x=results, file="model_ncut2_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_ncut2_error.pdf", plot=sim.model.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_ncut2_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_ncut2_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

sim.model.error=fsQCA.anchor.samford(raw.data=raw.data, depvar=c("raplib"), indepvars=c("manufac", "devalu", "hyperinf", "grostron", "groweak", "switcher", "execunco"), n.cut=3, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.65, given.incl.cut1=0.65, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.error$results
write.table(x=results, file="model_ncut3_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_ncut3_error.pdf", plot=sim.model.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_ncut3_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_ncut3_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

sim.model.error=fsQCA.anchor.samford(raw.data=raw.data, depvar=c("raplib"), indepvars=c("manufac", "devalu", "hyperinf", "grostron", "groweak", "switcher", "execunco"), n.cut=4, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.65, given.incl.cut1=0.65, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.error$results
write.table(x=results, file="model_ncut4_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_ncut4_error.pdf", plot=sim.model.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_ncut4_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_ncut4_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

sim.model.error=fsQCA.anchor.samford(raw.data=raw.data, depvar=c("raplib"), indepvars=c("manufac", "devalu", "hyperinf", "grostron", "groweak", "switcher", "execunco"), n.cut=5, reps=100, explain=1, error.reps=100, anchors=c(2), user.thresholds=user.thresholds, given.incl.cut0=0.65, given.incl.cut1=0.65, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.error$results
write.table(x=results, file="model_ncut5_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_ncut5_error.pdf", plot=sim.model.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_ncut5_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_ncut5_error_measure.pdf", plot=sensitivity.plot, scale=1.25)


### Simulate Model with Random Variable
sim.model.random=fsQCA.random(dataset=repdata.cal, depvar=c("raplib"), indepvars=c("manufac", "devalu", "hyperinf", "grostron", "groweak", "switcher", "execunco"), min.freq=1, max.freq=5, reps=100, explain=1, given.incl.cut0=0.65, given.incl.cut1=0.65, random.reps=100, min=.5)
write.table(x=sim.model.random$results, file="model_sim_random.csv", sep=",", row.names=FALSE)
percent=c(sim.model.random$percent.pos, sim.model.random$percent.neg, sim.model.random$percent.all)
percent.label=c("Percent Positive", "Percent Negative", "Percent All")
percents=data.frame(percent.label, percent)
write.table(x=percents, file="model_sim_random_percents.csv", sep=",", row.names=FALSE)
ggsave(filename="model_sim_random.pdf", plot=sim.model.random$plot, scale=1.85)


############################################################
### Replication & Simulation for Koenig-Archibugi (2004) ###
############################################################

### Set working sub-directory
setwd(main)
setwd("./Koenig-Archibugi (2004)")

### Import the data set
data <- read.csv(file="ka_data.csv")

### Isolate dependent and independent variables used in Model #1
repdata.cal=data[,c(2,3,5:7)]

### Set working sub-directory
setwd("./Results")

### Simulate Model #1
sim.model.1=fsQCA.sim(dataset=repdata.cal, depvar=c("supra"), indepvars=c("pubid", "conf", "reg", "mat"), min.freq=1, max.freq=10, reps=3000, explain=1, given.incl.cut0=-1, given.incl.cut1=-1, min=0)
write.table(x=sim.model.1$results, file="model_1_sim.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_sim.pdf", plot=sim.model.1$plot, scale=1.85)

### Isolate dependent and independent variables used in Model #2
repdata.cal=data[,c(2,4:7)]

### Simulate Model #2
sim.model.2=fsQCA.sim(dataset=repdata.cal, depvar=c("supra"), indepvars=c("opid", "conf", "reg", "mat"), min.freq=1, max.freq=10, reps=3000, explain=1, given.incl.cut0=-1, given.incl.cut1=-1, min=0)
write.table(x=sim.model.2$results, file="model_2_sim.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_sim.pdf", plot=sim.model.2$plot, scale=1.85)

### Isolate dependent and independent variables used in Model #1
repdata.cal=data[,c(2,3,5:7)]

### Simulate Model #1 with error (0, 0.01, 0.025, 0.05)
sim.model.error=fsQCA.anchor.ka(raw.data=repdata.cal, depvar=c("supra"), indepvars=c("pubid", "conf", "reg", "mat"), n.cut=1, reps=100, explain=1, error.reps=100, given.incl.cut0=-1, given.incl.cut1=-1, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.error$results
write.table(x=results, file="model_1_ncut1_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut1_error.pdf", plot=sim.model.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_1_ncut1_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut1_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

sim.model.error=fsQCA.anchor.ka(raw.data=repdata.cal, depvar=c("supra"), indepvars=c("pubid", "conf", "reg", "mat"), n.cut=2, reps=100, explain=1, error.reps=100, given.incl.cut0=-1, given.incl.cut1=-1, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.error$results
write.table(x=results, file="model_1_ncut2_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut2_error.pdf", plot=sim.model.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_1_ncut2_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_ncut2_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Isolate dependent and independent variables used in Model #2
repdata.cal=data[,c(2,4:7)]

### Simulate Model #2 with error (0, 0.01, 0.025, 0.05)
sim.model.error=fsQCA.anchor.ka(raw.data=repdata.cal, depvar=c("supra"), indepvars=c("opid", "conf", "reg", "mat"), n.cut=1, reps=100, explain=1, error.reps=100, given.incl.cut0=-1, given.incl.cut1=-1, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.error$results
write.table(x=results, file="model_2_ncut1_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut1_error.pdf", plot=sim.model.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_2_ncut1_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut1_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

sim.model.error=fsQCA.anchor.ka(raw.data=repdata.cal, depvar=c("supra"), indepvars=c("opid", "conf", "reg", "mat"), n.cut=2, reps=100, explain=1, error.reps=100, given.incl.cut0=-1, given.incl.cut1=-1, min=0, devs=c(0, 0.01, 0.025, 0.05))
results=sim.model.error$results
write.table(x=results, file="model_2_ncut2_error.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut2_error.pdf", plot=sim.model.error$plot, scale=1.85)
results["interval"]=cut(x=results[,1], breaks=100, labels=FALSE)
measure.melted=ddply(.data=results, .variables=.(interval, error), summarize, density=length(unique(Configurations))/length(Configurations))
sensitivity.measure=dcast(data=measure.melted,formula=interval~error )
sensitivity.plot=ggplot(data=measure.melted, aes(x=interval, y=density, colour=))+geom_smooth(aes(colour=error), size=1.5, se=FALSE)+xlab(label="Minimum Sufficiency Inclusion Score Percentile")+ylab(label="Unique Results Index")
write.table(x=sensitivity.measure, file="model_2_ncut2_error_measure.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_ncut2_error_measure.pdf", plot=sensitivity.plot, scale=1.25)

### Isolate dependent and independent variables used in Model #1
repdata.cal=data[,c(2,3,5:7)]

### Simulate Model 1 with Random Variable
sim.model.random=fsQCA.random(dataset=repdata.cal, depvar=c("supra"), indepvars=c("pubid", "conf", "reg", "mat"), min.freq=1, max.freq=2, reps=100, explain=1, given.incl.cut0=-1, given.incl.cut1=-1, random.reps=100, min=0)
write.table(x=sim.model.random$results, file="model_1_sim_random.csv", sep=",", row.names=FALSE)
percent=c(sim.model.random$percent.pos, sim.model.random$percent.neg, sim.model.random$percent.all)
percent.label=c("Percent Positive", "Percent Negative", "Percent All")
percents=data.frame(percent.label, percent)
write.table(x=percents, file="model_1_sim_random_percents.csv", sep=",", row.names=FALSE)
ggsave(filename="model_1_sim_random.pdf", plot=sim.model.random$plot, scale=1.85)

### Isolate dependent and independent variables used in Model #2
repdata.cal=data[,c(2,4:7)]

### Simulate Model 2 with Random Variable
sim.model.random=fsQCA.random(dataset=repdata.cal, depvar=c("supra"), indepvars=c("opid", "conf", "reg", "mat"), min.freq=1, max.freq=2, reps=100, explain=1, given.incl.cut0=-1, given.incl.cut1=-1, random.reps=100, min=0)
write.table(x=sim.model.random$results, file="model_2_sim_random.csv", sep=",", row.names=FALSE)
percent=c(sim.model.random$percent.pos, sim.model.random$percent.neg, sim.model.random$percent.all)
percent.label=c("Percent Positive", "Percent Negative", "Percent All")
percents=data.frame(percent.label, percent)
write.table(x=percents, file="model_2_sim_random_percents.csv", sep=",", row.names=FALSE)
ggsave(filename="model_2_sim_random.pdf", plot=sim.model.random$plot, scale=1.85)
